function putative_matches = getMatches(feat1, feat2, image1, image2, mws, sws, thresh)
% will return a matrix of 4 element vectos containing the matching vectors
    % feat* are the features or 'corners'. Must be in a nx2 = [x(n) y(n)]
    % image* are the images where featm = imagem and featn = imagen
    % mws is the 'matching window' size
    % sws is the search window size (pixels of the image to be considered in the
    % matching!)
    % thresh is a threshold for the simmilarity

%feat1=x_crn; feat2=xp_crn;mws= 9;sws= 120;thresh= 0.8 ;

%% First get the matches TO features in 'image1' WITH 'image2' and viceversa
%realtive mask
sz_img1 = size(image1);
sz_img2 = size(image2);
rel_mask = [(-mws/2 + .5) (-mws/2 + .5);(mws/2 + .5) (mws/2 + .5)]; %realtive positions given the compare mask window
matches_1_2 = getPotMatches(image1,image2,feat1,feat2,sws,thresh,sz_img1,sz_img2,rel_mask);%gets the cell
matches_2_1 = getPotMatches(image2,image1,feat2,feat1,sws,thresh,sz_img2,sz_img1,rel_mask);%gets the cell

%matches_1_2 will so contain the possible matches of corners in image1 with
%image 2. Similarly matches_2_1 will have matches of corners in image2 with
%image 1. Both are arranged in cell structures with matches_x_y{i} = Nx2
%where i is the index of the corner in image x and N is the number of
%possible matches with image 2. First row is the index of the specific
%corner, row 2 is the value of the correlation between points featx(i,:)
%and featy(row 1,:)

%% Second we get the best matches of each feature
for i = 1:length(feat1)
    if isempty(matches_1_2{i}) %if the current feature has no possible matches
        bm1(i) = length(feat2)+1; %an out of bounds index
    else
        [~,idx] = max( matches_1_2{i}(:,2) );%find maximum value of correlation
        idx = matches_1_2{i}(idx,1);%get index of the correlatad corner
        bm1(i) = idx;% index of the greatest match of ith feature of image one is stored in bm1(i)
    end
end
for i = 1:length(feat2)
    if isempty(matches_2_1{i}) %if the current feature has no possible matches
        bm2(i) = -1; 
    else
        [~,idx] = max( matches_2_1{i}(:,2) );%find maximum value of correlation
        idx = matches_2_1{i}(idx,1);%get index of the correlatad corner
        bm2(i) = idx;% index of the greatest match of ith feature of image 2 is stored in bm2(i)
    end
end
bm2(length(feat2)+1) = -1;
%now we have arrays bm1 and bm2. They hold the indexes of the best matches
%of each image. Now we find the matches that point to each other.
matchesIdx = find( bm2(bm1) == 1:length(feat1));
putative_matches = [matchesIdx' bm1(matchesIdx)'];


function mtchs = getPotMatches(image1,image2,feat1,feat2,sws,thresh,sz_img1,sz_img2,rel_mask)
mtchs = {};
for i = 1:length(feat1) %for each corner in image 1  
    this_feat = feat1(i,:); %get i-th corner in the (x,y) coord format
    search_area = [this_feat - sws/2; this_feat + sws/2];%the search area of other features
    %indices of features of image 2 inside the search area (reduces calculations):
    % TODO: make the search area have the same ratio of the image!
    idx_feats2inarea = find(feat2(:,1) > search_area(1,1) & ... x coord more than min x
                            feat2(:,1) < search_area(2,1) & ...
                            feat2(:,2) > search_area(1,2) & ...
                            feat2(:,2) < search_area(2,2) );
    msk = rel_mask + [this_feat;this_feat];%current coordinates according to rel_mask -> mws
    msk = [max(msk(1,1),1) max(msk(1,2),1);min(msk(2,1),sz_img1(2)) min(msk(2,2),sz_img1(1))]; %boundary handling
    pix_thisfeat_im1 = image1(msk(1,2):msk(2,2),msk(1,1):msk(2,1));
    pix_thisfeat_im1 = (pix_thisfeat_im1 - mean(pix_thisfeat_im1(:)));
    pix_thisfeat_im1 = pix_thisfeat_im1 / (norm(pix_thisfeat_im1(:))+eps);
    matches = idx_feats2inarea;
    for j = 1:length(idx_feats2inarea) %for each feature of image 2 inside the search area!
        idx = idx_feats2inarea(j);%current index of features from image2
        crn_thisfeat_im2 = feat2(idx,:);%current corner coords of image2
        msk = rel_mask + [crn_thisfeat_im2;crn_thisfeat_im2];%current coordinates according to rel_mask -> mws
        msk = [max(msk(1,1),1) max(msk(1,2),1); min(msk(2,1),sz_img2(2)) min(msk(2,2),sz_img2(1))]; %boundary handling
        pix_thisfeat_im2 = image2(msk(1,2):msk(2,2),msk(1,1):msk(2,1));%pixel mask of image2
        pix_thisfeat_im2 = (pix_thisfeat_im2 - mean(pix_thisfeat_im2(:)));
        pix_thisfeat_im2 = pix_thisfeat_im2 / (norm(pix_thisfeat_im2(:))+eps);
        if size(pix_thisfeat_im1) == size(pix_thisfeat_im2)
            corrl = dot(pix_thisfeat_im1(:),pix_thisfeat_im2(:));
        else
            corrl = -inf;
        end
        matches(j,1) =  idx;%so, the coords of the possible match will be in feat2(idx,:)
        matches(j,2) =  corrl; %will save indices and correlations.
    end
    if matches %if the variable matches is not empty
        %will save the potential matches of the ith feature of image 1 with
        %image 2, it saves them in a cell of size length(feat1)
        mtchs{i} = matches(matches(:,2) > thresh,:);
    else
        mtchs{i} = [];
    end
end







